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Abstract 

Background: Arsenic is an environmental pollutant, potent human toxicant, and oxidative stress agent with a 
multiplicity of health effects associated with both acute and chronic exposures. A semi-mechanistic cellular-level 
toxicokinetic (TK) model was developed in order to describe the uptake, biotransformation and clearance of 
arsenical species in human hepatocytes. Notable features of this model are the incorporation of arsenic-glutathione 
complex formation and a "switch-like" formulation to describe the antioxidant response of hepatocytes to arsenic 
exposure. 

Results: The cellular-level TK model applies mass action kinetics in order to predict the concentrations of trivalent 
and pentavalent arsenicals in hepatocytes. The model simulates uptake of arsenite OAs'") via aquaporin isozymes 9 
(AQP9s), glutathione (GSH) conjugation, methylation by arsenic methyltransferase (AS3MT), efflux through multidrug 
resistant proteins (MRPs) and the induced antioxidant response via thioredoxin reductase (TR) activity. The model 
was parameterized by optimization of model estimates for arsenite OAs" 1 ), monomethylated (MMA) and 
dimethylated (DMA) arsenicals concentrations with time-course experimental data in human hepatocytes for a time 
span of 48 hours, and dose-response data at 24 hours for a range of arsenite concentrations from 0.1 to 10 uM. 
Global sensitivity analysis of the model showed that at low doses the transport parameters had a dominant role, 
whereas at higher doses the biotransformation parameters were the most significant. A parametric comparison of 
the TK model with an analogous model developed for rat hepatocytes from the literature demonstrated that the 
biotransformation of arsenite (e.g. GSH conjugation) has a large role in explaining the variation in methylation 
between rats and humans. 

Conclusions: The cellular-level TK model captures the temporal modes of arsenical accumulation in human 
hepatocytes. It highlighted the key biological processes that influence arsenic metabolism by explicitly modelling 
the metabolic network of GSH-adducts formation. The parametric comparison with the TK model developed for 
rats suggests that the variability in GSH conjugation could have an important role in inter-species variability of 
arsenical methylation. The TK model can be incorporated into larger-scale physiologically based toxicokinetic (PBTK) 
models of arsenic for improving the estimates of PBTK model parameters. 
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Background 

Arsenic is a naturally occurring metalloid, abundant in 
the earth's crust and a component of more than 245 
minerals [1]. Exposure to arsenic has been associated 
with cancers of the liver, bladder, skin and lung [2,3]. 
Epidemiological studies in Taiwan, Bangladesh and India 
have reported adverse health effects associated with 
chronic arsenic exposure including; chronic obstructive 
pulmonary disease, non-cirrhotic portal fibrosis, hyper- 
tension and ischeamic heart disease [4], The risk of 
developing serious diseases from chronic exposure to 
inorganic arsenic in drinking water prompted the US 
Environmental Protection Agency (EPA) to lower the 
maximum contamination level (MCL) for arsenic in 
drinking water to 10 ppb [5]. 

There are two biologically important arsenic valence 
states: arsenite (As(OH) 3 , iAs m ) and arsenate (AsO(OH) 3 , 
iAs v ). Inorganic arsenic in water is largely in the form of 
arsenate; it is negatively charged at physiological pH and 
slowly taken up by cells [6] . Arsenate is rapidly converted 
to arsenite in vivo [7] which is taken up by cells much 
more quickly than arsenate [8] . Methylation of arsenicals 
facilitates their excretion from the cell and therefore was 
long considered a detoxification process, but recent evi- 
dence indicates that monomethylated (MMA) and 
dimethylated (DMA) arsenicals have many toxic effects 
including increased oxidative stress [9], chromosomal 
aberrations (CA), and oxidative DNA damage [10-12]. In 
hepatocytes, trivalent monomethylated arsenicals 
(MMA 111 ) inhibit the activity of thioredoxin reductase 
(TR), which is a critical antioxidant enzyme controlling 
the cellular redox balance [13,14]. 

Uptake and efflux of arsenicals occur primarily through 
transporter proteins. Uptake of iAs m in hepatocytes and 
efflux of MMA 111 to blood take place through aquaporin 
isozymes 9 (AQP9), a family of membrane-spanning pro- 
teins that facilitate movement of solutes down their con- 
centration gradient. AQP9 channels are expressed at high 
concentrations in liver cells and have been shown to trans- 
port iAs 111 when expressed both in Saccharomyces cerevi- 
siae (yeast) and in Xenopus oocytes [15-17]. Another class 
of transmembrane proteins facilitating the transport of 
iAs 111 and MMA 111 across the cellular membrane of hepa- 
tocytes is glucose transporters and especially GLUT2 
which is highly expressed in the liver [18,19]. Glutathione 
conjugated arsenicals are exported to the extracellular 
space via multidrug resistant proteins (MRPs) and multi- 
drug resistant P-glycoproteins (PGPs)which are ATP- 
binding cassette (ABC) transporters that export solutes 
against their concentration gradient [19-21]. 

Methylation of inorganic arsenic takes place primarily in 
the liver and specifically in hepatocytes via enzymatic cata- 
lysis by arsenic methyltransferase (AS3MT), previously 
known as Cytl9, producing both mono- and dimethylated 



arsenicals [22-25]. Two biochemical pathways have been 
proposed for methylation of arsenates with a key differ- 
ence in the substrate for AS3MT methylation: (a) a classi- 
cal process of alternating steps of reduction and oxidative 
methylation where iAs and MMA are the substrates 
and the methylation can happen in the presence or 
absence of GSH [26] and (b) a process of GSH conjugation 
and reductive methylation where thiol-bound arsenicals 
(arsenic triglutathione - ATG, monomethylarsenic diglu- 
tathione - MADG) are the substrates [27,28]. GSH has a 
stimulatory role in both methylation pathways either as a 
reductant or in direct conjugation with arsenicals [29]. 

Arsenic activates the redox sensitive transcription fac- 
tor Nuclear Factor -E2- related factor 2 (Nrf2) causing 
its increased nuclear translocation and Nrf2 binding to 
the Antioxidant Response Element (ARE) [30,31]. 
Arsenic activates Nrf2 in a different manner when com- 
pared to other compounds such as sulforaphane (SF) 
and tert-butylhydroquinone (tBHQ) enhancing the inter- 
action of specific subunits of the E3 ubiquitin ligase 
[32]. It has been suggested that hepatocytes exhibit an 
"on-switch" antioxidant behavior when exposed to 
increasing arsenic doses [33], possibly a result of this 
Nrf2 activation. 

In this study, a cellular-level semi-mechanistic TK 
model was developed for predicting intra-cellular con- 
centrations of different arsenicals (trivalent and pentava- 
lent) in hepatocytes. Currently, the only published 
cellular-level TK model for the uptake, biotransforma- 
tion and efflux of arsenicals is the Easterling et al. [34] 
model schematically shown in Figure 1; they demon- 
strated the relative importance of transport processes 
affecting the accumulation of arsenicals in rat hepato- 
cytes [34]. However, a TK model for humans is needed 
since the inherent variation of arsenic metabolic capa- 
city across various organisms complicates the extrapola- 
tion of rat TK model to humans. Therefore, the 
mathematical model presented here was developed for 
humans and parameterized based on data from human 
hepatocytes. 

Methods 

The human TK model applies mass action kinetics in 
order to predict the concentrations of trivalent and 
pentavalent arsenicals including arsenite (iAs ), 
monomethylated (MMA), and dimethylated arsenicals 
(DMA) in human hepatocytes. This TK model takes 
into account processes such as influx, efflux, methyla- 
tion, oxidation and glutathione conjugation of arseni- 
cals. Moreover, it accounts for induced cellular 
antioxidant response due to arsenic exposure through 
a "switch-like" mechanism that alters the model 
response above a specific threshold concentration 
[35]. 
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Figure 1 Schematic depiction of the cellular level TK model for 
rat hepatocytes. The solid lines that cross the cellular membrane 
(oval) represent transport processes, while the solid lines within the 
oval represent biotransformation. The dashed line represents the 
inhibitory effect of iAs'" on the second methylation reaction (MMA 
to DMA) (Source: Easterling et al. [34]) 



This model has been compared with the Easterling 
et al. model [34], in terms of the ability to fit to data of 
arsenic retention and methylation in human hepato- 
cytes. This comparison aims to highlight the advantages 
of developing biologically relevant TK models based on 
data acquired from human cells. Further comparison of 
these two models in terms of their estimated parameter 
values aim to study the major intracellular kinetic pro- 
cesses that contribute to the differences in metabolism 
between humans and rats. 

Model Development 

The semi-mechanistic TK model describes arsenic trans- 
port across the cellular membrane and arsenic metabo- 
lism in hepatocytes according to the metabolic reaction 
cascade proposed by Hayakawa et al. [27] (Figure 2). 
Figure 3 presents a schematic depiction of the constitu- 
ents of the TK model that are also explained in Table 1. 
These constituents include chemical species and 
enzymes, and the interactions among them. Fundamental 
assumptions made in the formulation of the TK model 
are: 
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Figure 2 A new metabolic pathway of inorganic arsenic 
biotransformation via arsenic-GSH complexes formation This 
pathway includes two separate branches of arsenic 
biotransformation: MADG -> MMA 1 " -> MMA V and MADG -> DMAG 
-> DMA 1 " -> DMA V (Source: Hayakaya ef al [27]). 
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Figure 3 Schematic depiction of the cellular level TK model for 
human hepatocytes showing the major components. Green 
squares represent extracellular amounts of arsenicals, while red 
squares represent intracellular amounts. Ovals represent activities of 
proteins (AQP9, TR, MT1, MT2, MRPa and MRm) and GSH (GS-Pm, 
GS-Pd). Arrows and hammerheads indicate activation and inhibition 
respectively. 
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Table 1 Optimized parameter values of the TK model along with the corresponding process they describe 



Notation 


Parameter 


Value 


Units 


Description 


1 


K ATG int 


0.25 


1/min 


Rate of production of ATG 


2 


L, 

■wlADG int 


90 


1/min 


Rate of production of MADG 


3 


I 

K DMAG 


0.01 22 


1/min 


Rate of production of DMAG 


4 




33.254 


1/min 


Rate of oxidation of DMA 1 " 


5 


L 

^oxm 


0.2375 


1/min 


Rate of oxidation of MMA 1 " 


6 


i III 

K iAs int 


320 


1/min 


Rate of production of iAs 1 " 


7 


i HI 
K MMA int 


1 200 


1/min 


Rate of production of MMA I!I 


8 




8 




Hill coefficient of inhibition of MMA 1 " production 


9 


KU .; 


1 2.94 


|jM 


Dissociation constant of inhibition of MMA 111 


10 


f 

T GSHm 


0.992 




Coefficient of inhibition of MADG hydrolysis 


1 1 


Kgin 


704.96 


1/min 


Rate of GSH production increase 


12 


1, III 
K DMA int 


0.8472 


1/min 


Rate of production of DMA 1 " 


13 


r 

r GSHd 


r\ r\r\o o 

0.9988 




Coefficient of inhibition of DMAG hydrolysis 


1 4 


Vmax t2 


1 .237 


pmol/min 


Maximal rate of MADG efflux 


15 


Km t2 


1 9.47 


l-iiVI 


Half saturation constant of MADG efflux 


16 


km in 


512.27 


1/min 


Rate of upregulation of MRP 


1 7 


1, 

k ss 


4.26 


1/min 


Steady state rate of efflux of iAs 1 " 


18 


T e 


1 0 


m i n 


Time constant of AQP9 inactivation 


19 


l , 


4.2 


1/min 


Initial rate of efflux of iAs 1 " 


20 


1, 

Kinf 


1 .61 3 


1/min 


1— .£1, ,., ;a^IH 

Intlux ot iAs 


21 


t, 

K MMAext 


0.0006 


1/min 


Rate of efflux of MMA 


22 


r 

fm 


0.2 




Coefficient of efflux of MMA 


23 


Vmax tq 


0.35 


pmol/min 


Maximal rate of ATG efflux 


24 


Km tl 


2.3 


uM 


Half saturation constant of ATG efflux 


25 


L 

K DMAext 


0.002 


1/min 


Kate ot ettlux ot DMA 


26 


r 

td 


3 




Coefficient of efflux of DMA 


27 


Vmax m 


50 


pmol/min 


Maximal rate of ATG methylation 


28 


Km m 


9.32 


uM 


Half saturation constant of ATG methylation 


29 


ni 


1 .22 




Hill coefficient of ATG methylation 


30 


Kd 


0.31 5 


uM 


Dissociation constant of ATG methylation 


31 


Kim 


1 .53 


(jM 


Inhibition constant of ATG methylation 


32 


Kid 


1 


uM 


Inhibition constant of MADG methylation 


33 


VmaXd 


8 


pmol/min 


Maximal rate of MADG methylation 


34 


Km d 


0.034 


(jM 


Half saturation constant of MADG methylation 


35 


n 2 


1 .83 




Hill coefficient of MADG methylation 


36 


Kd i 


2 33 


(jM 


lii c c r~\r i a t i i"^ n rnnctant r~ir Iv /I A 1 mothu atinn 
Ul SSULId UUi 1 LUIIbldllL Ul IVIALyVj 1 1 Itrll ly Id UUI 1 


37 


ka in 


1.64 


1/min 


Rate of AS3MT efficiency increase 


38 


f A 


50 




Coefficient of MADG methylation inactivation 


39 




0.99 


1/min 


Rate of TR signaling 


40 


^TRinact 


987.13 


1/min 


Rate of TR inactivation 


41 


TRc 


15 




Constant of TR inactivation 


42 


N 


1.75 




Hill coefficient of TR inactivation 


43 


TR 0 


9.27 


1/min 


Steady state value of TR activity 
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1) Arsenite influx across the cellular membrane 
pores is governed by their electrochemical potential, 
and can be described through an ion channel con- 
ductance-based formulation. 

2) The oxidative stress mediated response of hepato- 
cytes to arsenite exposure exhibits a "switch-like" 
behavior, and the upregulation of enzyme activities 
can be described through an approximate step func- 
tion at a threshold concentration. 

3) The methylation reactions are influenced by coop- 
erativity phenomena as well as substrate inhibition, 
and can be described through a hybrid approach of 
Hill and Michaelis-Menten kinetics. 

4) The GSH-bound hydrolysis and clearance of 
methylated arsenicals exhibit a threshold-dependent 
behavior, and can be described using a sigmoidal 
function. 

5) Concentrations of arsenicals are uniform within 
the hepatocytes as well as the extracellular medium 

6) All hepatocytes in the system have identical prop- 
erties, are uniformly distributed in the medium, and 
are exposed to the same extracellular concentrations 
of arsenicals. 

Uptake of arsenite by hepatocytes via AQP9s [36] is 
governed by their electrochemical potential across the 
cellular membrane (for simplicity we refer to AQP9s as 
being the ensemble of the activity of both AQP9 and 
GLUT2 channels). The conductance-based formulation 
for ion channels proposed in the Hodgkin-Huxley 
model [37] is used here to describe the regulation of 
arsenite flux by AQP9s (Equations 1-2). Specifically, the 
inactivation of AQP9 subunit gates during iAs 111 influx 
is described by Equation lb, which expresses the 
increased probability of these gates closing as more 
transport across the gates occurs. 



d(iAs'") i 
dt 



= k inf x (iAs m ) m - AQP9 x (iAs 111 ) 



(la) 



k ATG . x(iAs m ) int + k m x (ATG); 

int iAs im 



family of enzymes; is the rate constant for 

hydrolysis of ATG (reciprocal to conjugation); kss is the 
steady state rate constant for efflux of arsenite that is 
attained at long time periods; k 0 is the rate constant 
describing the basal activity of AQP9; and T e is the time 
constant governing the regulation of AQP9 gates. 

Thioredoxin (Trx) Reductase (TR) is the enzyme that 
catalyzes the reduction of Trx. Trx is a critical antioxi- 
dant protein and an important reductant in the methyla- 
tion of arsenic by AS3MT [38]. The inactivation of TR 
by MMA m leads to signals that account for two differ- 
ent phenomena: the induction of GSH and ABC trans- 
porters via a redox sensitive activation of the cellular 
antioxidant response Nrf2 nuclear receptor pathway 
[14,39] and the decreased methylation capacity of 
AS3MT. In this study the main focus is on MRPs as an 
efflux mechanism of arsenicals since it has been 
reported that they are regulated by Nrf2 [40,41]. The 
inactivation is modeled using principles of indirect 
response model theory [42,43] via the threshold-depen- 
dent parameter: s = [iAs_J ja j L:i threshold Parameter s 

threshold 

depends on the initial exposure concentration of 
arsenite, (iAs nl ) init and a threshold concentration. The 
value of S is zero when arsenite doses are below the 
threshold concentration, and gradually increases with 
greater arsenite doses. The following equations describe 
this reaction cascade: 



dTR 

xk 



k TR x(TR 0 -TR) 



TR, 



TR;„ 



xH 



TR 



111 iN 



H 



[MMA 

TR ~ (IC TR ) N + [MMA m ] N 



dGSH . 1 

= k„ x k„ xGSH 

dt gin TR g deg 



(3a) 



(3b) 



(4a) 



AQP9 = [k 0 +(k ss -k 0 )x(l-e ^)] 



(lb) 



g deg TR, 



(4b) 



d(iAs"') e 
dt 



= AQP9x (iAs ) int - k inf x (iAs m ) e 



(2) 



Where, k inf represents the mass transfer coefficient for 

influx of arsenite in hepatocytes; ^ATG int represents the 

rate constant for arsenite conjugation with GSH to form 
ATG catalyzed by the Glutathione S-Transferase (GST) 



dMRP , 1 

= k m x k m xMRP 

dt ^ TR m deg 



deg TR r 



(5a) 



(5b) 
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dAS3MT 

dt 



xTR-k 3 x AS3MT 

a dec 



d dec 



X TR r 



(6a) 



(6b) 



where, k TR is the first-order rate constant controlling 
the activity of TR; TR 0 is the baseline activity value of 

TR; k TR is the first-order rate constant for TR 

inact 

inactivation; TR C is a dimensionless inactivation con- 
stant; N is the Hill coefficient for enzyme inactivation 

from MMA m ; k„ k m . and k a . are the rate con- 

6 in in in 

stants governing of the activities of GSH, MRP and 

AS3MT, kg deg , ^ m deg an(jl ^ a dec are t ^ ie corres P onc '- 
ing first-order decay constants. 

Methylation reactions of arsenic in the liver have been 
modeled in published cellular-level and whole-body 
PBTK models with classical Michaelis-Menten kinetics 
[34,44-46]. Alternatives include cooperativity models such 
as the classical Hill-type formulation and the more 
mechanistic Monod-Wyman-Changeux (MWC) model 
[47]. In preliminary analyses, these formulations were 
unable to explain the time course patterns of arsenic 
retention and methylation in human hepatocytes (results 
not shown). Therefore, an alternative, non-classical formu- 
lation was used in the TK model. In this model, the 
AS3MT is assumed to exhibit cooperativity and the 
methylation reaction rate is assumed to exhibit hysteretic 
sigmoidal behavior as per Frieden [48]. In this formulation, 
the cooperativity is described by a Hill-type formulation 
for V max that is dependent on the total ATG present in 
the system. The formulation accounts for the constitutive 
influence of GSH in the methylation reaction cascade and 
its role in the increase of V max . Moreover this Hill-type 
formulation for V max accounts for a possible colocalization 
of thiol-containing proteins that interact with GSH (e.g. 
GSTP1), MRPs and AS3MT in hepatocytes. This colocali- 
zation would control not only the production and clear- 
ance of ATG but methylation activity as well [49,50]. 

Previous studies have suggested that exposure of 
human hepatocytes to elevated doses of iAs (0.4 - 
4 uM) markedly reduced the production of DMAs while 
at the same time increased the yields of MMAs [51]. 
Therefore, it is assumed here that the AS3MT inactiva- 
tion signal (Equation 6) affects only the second methyla- 
tion reaction rate (Equation 8a). 



MADG, 



MT 



[ATG]; 



1 + 



[ATG], 
Km m 



x(l + 



[ATG]; 



(7a) 



MT, 



([ATG] int + [ATG] 

ext J 

(KdJ n i+([ATG] int +[ATG] ext ) n i 

Vmax m 



(7b) 



DMAG, 



•■ MT, 



[MADG]; 



1 + 



[MADG]; 
Km A 



x(l + 



[ATG]; 



(8a) 



MT, 



AS3MT 



f A x tanh(S) + 1 

([ATG] int +[ATG] ext ) n 2 ^ Vmax 

(Kd 2 ) n 2+([ATG] int +[ATG] ext ) n 2 Km c 



(8b) 



where, MADG ui and DMAG ui are the rates of arsenic 
methylation for the first and second methylation reac- 
tions respectively. ^max m and V maXd are the maximal 
rates of the first and second methylation reactions 
respectively; ^m m and Km d are the half-saturation 
constants for the methylation reactions; K; m and Kid are 
the uncompetitive inhibition constants for the respective 
reactions ni, n 2 are the Hill coefficients; K d and K d2 
are the dissociation constants influencing the sigmoidal 
change in ^max m and V maXd ; and f A is a coefficient of 
the second methylation reaction inactivation. 

MADG hydrolysis reaction is modeled using a 
"switch-like" formulation. For doses below the threshold, 
a Hill-type formulation is used (Equation 9b). Above the 
threshold concentration, the rates of hydrolysis of 
MADG and DMAG to MMA 1 " and DMA 111 , respec- 
tively, are assumed to be attenuated due to oxidative 
stress-induced GSH upregulation (Equations 9-10). On 
the other hand, this non-linear behavior may result in 
an increase or a decrease of GSH depending on the con- 
centration of iAs and the duration of exposure. The 
non-linear sigmoid function tanh(S) is used here to 
describe this "switch-like" behavior; this formula has 
been previously used [52] in a neurocomputational 
model to describe the non-linear threshold-dependent 
behavior of neuronal firing rate. 



HD m = [l-f GSHm Xtanh(S)] 

x (l-H GSH + tanh(S)xH GSH ) 
xGS-P m x(MADG) int 



MMA 

GS-P m - lnt 



(9a) 



(9b) 



GSH 



Stamatelos et al. BMC Systems Biology 201 1, 5:16 
http://www.biomedcentral.eom/1 752-0509/5/1 6 



Page 7 of 15 



II 



GSH 



([ATG] int + [ATG] 

ext J 



(Kd 3 ) n 3+([ATG] int +[ATG] ext ) 1 



HD d = [l-f GSHd Xtanh(S)]xGS-P d 
x (DMAG) int 



(9c) 



(10a) 



d(MMA'") e 
dt 



[(f m -l)xtanh(S) + l] 



x k_ 

d(DMA) v 
dt~~ 

X knyA 



MMA x (MMA m )j 

ext 



- = [(f d -l)xtanh(S) + l] 
x (DMA v ) int 



(14) 



(15) 



GS-P, = 



DMA 



(10b) 



GSH 



where HD m and HD d are the rates of hydrolysis of 



MADG and DMAG respectively. * 



MMA 



is the reac- 



tion rate constant for MADG hydrolysis (MMA 111 pro- 
duction); fGSH is the coefficient of inhibition of 

m 

MADG hydrolysis; n 3 and K d3 are the Hill coefficient 
and dissociation constant, respectively, for the inhibition 



term; 



,111 



the reaction rate constant for 



DMAG hydrolysis (DMA 111 production); and ^GSH d is 

the coefficient of inhibition of DMAG hydrolysis. 

Efflux of GSH (or protein) -bound arsenic adducts 
(ATG, MADG) is assumed to take place via multidrug 
resistant proteins (MRPs), and is described by classical 
Michaelis-Menten kinetics. Since MADG is a substrate 
in the dimethylation reaction (Equation 8a), its efflux 
rate is assumed to be affected by MRP levels [51]. 



d(ATG) e 
dt 



MRP, x 



[ATG]; 



1 + 



[ATG], 

Km* 



Vn 



MRP, 



Km, 



(11a) 



(lib) 



d(DMA) 



in 



dt 



:[(f d -l)xtanh(S) + l] 



(16) 



xk DMAext X[DMA m )i 



where, ^max '^r 



and 



V„ 



are the 

Michaelis constants of the biophysical clearance of ATG 

and MADG, respectively; ^mma and k DMA are 

the rate constants of MMA and DMA clearance, respec- 
tively; f m and f d are the dimensionless coefficients of 
clearance of the respective processes affecting the maxi- 
mal efflux. 

The remaining biotransformation reactions include a 
series of methylation, glutathione conjugation and oxi- 
dation reactions [27] (Equations 17-23). 



d(ATG), 
dt~~ 



k ATG x(iAs m ) 



iAs ' 



x (ATG) int - MADG t 



d(ATG) e 
dt 



d(MADG) int 
dt 



MADG, 



+ kMADG. x(MMA m ) int 

int 

d(MADG) ext 



DMAG,,: - HD r 



dt 



(17) 



(18) 



d(MADG) ext =MRp x [MADG] 



dt 



Vmax t 

MRP m = MRP x ■ - 



1 | [MADG] int 

Km, 



Km, 



d(MMA ) f 
dt 



:[(f m -l)xtanh(S) + l] 



xk M MA . x (MMA )j 

ext 



(12a) 



(12b) 



(13) 



d(MMA m ) 
dt 

x (MMA 1,1 - 
d(MMA v )j 



HDm (kMADG jn , + k OX m) 

„, , d(MMA III ) ext 



dt 



x (MMA 111 ) 



dt 

d(MMA v ) e 
dt 



(19) 



(20) 
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d(DMAG) 

dt 



DMAG,,; + k 



DMAG 



(21) 



x(DMA m ) int -HD d 



d(DMA ) 
dt 



HD d - (k DMAG + k oxd ) 



(22) 



d(DMA v ) int _ 
dt 

d(DMA) x 



dt 



^oxd x (DMA ); 



(23) 



dt 



where kf^^Q is the rate constant of MADG pro- 
mt 

duction catalyzed by the Glutathione S-Transferase 
(GST) family of enzymes; k oxm is the rate constant 
for MMA m oxidation; koMAG is the rate constant for 
DMAG production catalyzed by the Glutathione 
S-Transferase (GST) family of enzymes and k oxl j is the 
rate constant for DMA oxidation. It has been sug- 
gested for this biotransformation pathway that trivalent 
arsenicals mostly bound to thiol-containing proteins are 
conjugated with GSH and methylated in the presence of 
arsenic methyltransferase (AS3MT) [28]. Therefore, the 

parameters k A TGint> ^MADG int and k DMAG indirectly 

represent the binding of trivalent arsenicals to thiol- 
containing proteins in the TK model. 

The TK model has been implemented in MATLAB; 
the system of ODEs comprising the TK model is 
solved numerically using the stiff solver odelSs. First, 
the model parameters corresponding to low doses (i.e. 
below the threshold) were estimated using time course 
in vitro measurements of arsenicals following expo- 
sures to 0.1 uM iAs 111 , from Styblo et al. [51]; it was 
assumed that at this dose, the hepatocytes exhibit no 
induced antioxidant response. Subsequently, model 
parameters corresponding to a wider range of doses 
(i.e. including both low-dose and high-dose behavior) 
were estimated using dose-response data (for doses 
ranging from 0.1 - 10 uM) reported by Drobna et al. 
[53]. The Drobna dataset includes measured concen- 
trations of iAs 111 , MMA, and DMA in primary cultured 
human hepatocytes after 24 h exposure to iAs (data 
for hepatocytes from 8 donors). For this case study, 
data on hepatocytes from one donor (white female, 
aged in the 60 s, Donor C) [53] were used, since this 
donor had similar characteristics to the human donor 
in the study by Styblo et al. [51]. The deterministic 



optimization function fmincon was used for parameter 
estimation in both cases. 

Sensitivity analysis 

Sensitivity analysis provides estimates of how variation 
of model's output can be apportioned to different 
sources of variation in model parameters. This quantity 
is given by the formula: 



var P [E(Y/P)] _ Dp 
var(Y) D 



(24) 



where Y denotes model output and P denotes the vec- 
tor of model parameters; D P and D illustrate the partial 
and total variance of the model output due to variation 
in model parameters according to assigned statistical 
distributions. 

The Fourier Amplitude Sensitivity Test (FAST) 
decomposes the total variance of model output (D) into 
terms of increasing dimensionality. FAST computes the 
Total Sensitivity Indices (TSI), which account not only 
for the variance due to individual parameters (Dj), but 
also estimate the variance due to interaction among 
parameters (D ij) D^, etc.) The total variance for n 
dimensions is given by 



D : 



Z D i+ZZ Di ' + - +Di - 



2,..,n 



(25) 



i=l i=i+l 



The model parameters were assumed to be normally 
distributed with a coefficient of variation up to 10%; for 
some parameters, the coefficient of variation was 1%. Ten 
thousand (10,000) samples were generated and the nor- 
mal distributions for all parameters were truncated at 1% 
and 99% (approximately three standard deviations from 
the mean value). Three model outputs were selected for 
the sensitivity analysis: Areas under the Curve (AUCs) of 
total retention of MMA, DMA and iAs 111 in human hepa- 
tocytes. The SIMLAB modeling platform [54] was used 
to perform the global sensitivity analysis. 

Comparison with the TK modeling formulation for rat 
hepatocytes 

The TK model presented here was also compared with 
results from a published model for rat hepatocytes, in 
order to assess the inter-species differences and the fea- 
sibility of direct, cross-species extrapolation. Specifically, 
the comparison focused on major intracellular processes 
that influence the different metabolizing rates between 
these two species. First, the TK modeling formulation of 
Easterling et al. [34] was parameterized using data from 
Styblo et al. [51]. Although, a direct comparison is not 
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possible because of major differences in model struc- 
tures, a subset of parameters was selected for compari- 
son that describe three major biochemical processes 
which account for similar cellular phenomena in both 
models. These processes include transport of arsenite 
across cellular membrane, methylation of arsenic, and 
biotransformation of AS3MT substrate (iAs in the rat 
model, and ATG in human TK model). Specific para- 
meters used in comparison include (a) normalized activ- 
ity of AQP9 (NP), which is defined as the ratio of influx 
and efflux of arsenite, (b) normalized activity of AS3MT 
(NM), which is defined as the ratio of the corresponding 
reaction parameters, and (c) bioavailability of AS3MT 
substrate (BMRS), which is defined as the ratio of rates 
of hydrolysis and conjugation of iAs for the human 
model, and the dissociation constant for the protein 
binding of arsenic in the rat model. These parameters 
are specified by Equation 26, and are described in 
Table 2. Parameters governing the efflux of arsenic were 
not compared because there is no direct correspondence 
between the modeling formulations. 

In order to facilitate direct comparison, approximate 
volumes of the cell cultures used in different experi- 
ments (per well) were estimated. Hepatocytes were 
assumed to have a spherical shape with 25 um diameter 
for both humans and rats [55]. Cellular volumes (Vole) 
for human and rat were estimated to be 1.6 and 0.8 uL 
respectively, based on the number of cells used in 
human hepatocytes experiments (2*10 ) and rat hepato- 
cytes experiments (10 5 ) [34,51]. 



Table 2 Comparison of selected parameter values 
between two hepatocyte-level TK modeling formulations 



Human TK : 



NP: 



Mnf 



NM (1 st ) 



NM (2 nd ) 



Vn 



Km m X Vok 
Vmax^ 
Km d X Vok 



(26) 



ATG 

BMRS = 525- 



Rat TK: 



NP = 



K t x Volcxp 1 
NM (1 st ) = k 3 

NM (2 nd ) Vmax 



K m xVolc 



BMRS = ^ 
ki 



Kinetic Process 


Human-Hepatocyte 
TK model 


Rat-Hepatocyte 
TK model 


NP - iAs 1 " 


0.38 


0.38 


NM - 1 st Reaction (1/min) 


3.35 


0.02 


NM - 2 nd Reaction (1/min) 


147.06 


0.17 


Biotransformation - 1 st MRS 


0.0008 


33.33 



The two formulations are developed based on data of arsenic uptake and 
biotransformation in human and rat hepatocytes (NP: Normalized AQP9 
activity, NM: Normalized AS3MT activity, MRS: Methylation Reaction Substrate). 



Results and Discussion 

The semi-mechanistic TK model was parameterized 
using the fmincon function in MATLAB and time 
course data of arsenicals in human hepatocytes from 
Styblo et al. [51]. The parameters are shown in Table 1. 
The model was able to capture the three distinct modes 
of the time course patterns corresponding to experimen- 
tal data (shown in Figure 4, Row 1). In contrast, the 
only currently existing cellular level TK model for 
arsenic, from Easterling et al. [34], was parameterized 
using the same optimization technique and data, but the 
model was not able to adequately capture these modes 
(shown in Figure 4, Row 2). 

The time-course estimates from the TK model show 
that initially (within first minutes of exposure) the rate of 
influx of AQP9s is substantially higher than the metabo- 
lism, thus leading to a fast accumulation of arsenite 
inside the cells. Then, the influx is reduced, and metabo- 
lism increases, thus leading to a slow decrease in arsenite 
levels (till 8 to 9 hours). During this period, MMA pro- 
duction appears to be the dominant process, as shown by 
higher levels of MMA compared to DMA, attributable to 
the high rate of MADG hydrolysis. Subsequently, the 
arsenite concentrations decrease at a faster rate, the sec- 
ond methylation reaction becomes more significant, and 
MADG hydrolysis is inhibited (Equation 9). 

Figure 5 shows the dose-response profiles estimated 
by the TK model parameterized using Drobna et al. 
data [53]. The model explained the dose-response pro- 
files in the data, and captured the significant decrease in 
DMA amounts at higher arsenite doses. Based on the 
sensitivity testing the threshold concentration value of 
0.1 uM was able to adequately explain the arsenicals 
retention and metabolism, as shown in Figure 6. 
Threshold values above 0.1 uM overestimate the con- 
centration of DMA in hepatocytes by one order of mag- 
nitude in the low dose region. On the other hand, 
threshold values below 0.1 uM (e.g. 0.01 uM) underesti- 
mate the DMA concentration substantially. 

Results of the sensitivity analysis showed that the rela- 
tive contribution of variance of individual TK model 
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Figure 4 Predicted time course profiles of arsenicals in human hepatocytes using two TK model formulations. The top row shows 
estimates from the TK model presented here, while the bottom row shows estimates from the TK model formulated using the Easterling et al. 
[34] approach. The left column shows intracellular levels of iAs'", total MMA, and total DMA, while the right column shows corresponding 
extracellular levels. Experimental data are from Styblo ef al. [51] (exposures to 0.1 uM arsenic for 2 days). 



parameters varied significantly across different doses of 
arsenite. As shown in Figure 7A, at low doses (0.1 uM), 
the transport parameter k 0 (Parameter 19) contributes 
most to the variance in intracellular MMA levels. This 
agrees with Easterling et al. [34], who reported that the 
transport parameters are the most significant in relation 
to intracellular concentration of arsenicals. On the other 
hand, at higher doses, the parameters related to intracel- 
lular biotransformation of MMA are the most influen- 
tial. For 1 uM dose of iAs the most significant 
parameters are k oxm (Parameter 5), ^GSH m (Parameter 
10), ^MMA 111 . (Parameter 7) and K m t ^ (Parameter 15). 
The first three parameters (Parameters 5, 10, and 7) 
directly influence oxidation and glutathione conjugation 
reactions involving MMA , whereas the Michaelis 



constant (Parameter 15) controls the activity of MRPs 
that efflux MADG from the cells. At high doses, 
induced antioxidant response of hepatocytes to arsenic 
leads to increased production of GSH and MRPs in the 
cells, leading to higher production of MADG, which can 
be readily effluxed via membrane-associated proteins. 
On the other hand, when the oxidation reaction is 
dominant (k oxm ; Parameter 5), it results in higher pro- 
duction of MMA V , which becomes accumulated in the 
cells, leading to overall increase in intracellular MMA 
levels. 

Figure 7B shows sensitivity analysis results for intra- 
cellular DMA levels. Similar to results from sensitivity 
analysis of MMA levels (Figure 7A), the transport para- 
meter k 0 (Parameter 19) is the most influential 
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Figure 5 Predicted dose-response profiles in human hepatocytes using the cellular level human TK model The left panel shows total 
amounts (in hepatocytes and the medium) of iAs'", MMA and DMA. The right panel shows corresponding intracellular levels. Experimental data 
are from Drobna ef al. [53] for hepatocytes from a 63 year old white female. 



parameter at low doses (0.1 (iM). For 1 uM dose, DMA 
production is significantly influenced primarily by ^GSH d 
(Parameter 13), which affects the rate of DMAG hydro- 
lysis. Furthermore, the oxidation reaction is not impor- 
tant in the case of DMA levels, because DMA transport 
across the cellular membrane is much faster compared 
to MMA transport [56]. This is also corroborated by the 
relative values of corresponding transport parameters 
for MMA and DMA (f d >> f m> , as shown in Table 1), 
and by the findings of Styblo et al. [51]. 



Figure 7C shows sensitivity analysis results for intra- 
cellular iAs 111 levels. At low doses (0.1 uM), the TSIs of 
most parameters are close to 1, indicating very high 
contributions. This unusual finding can be attributed to 
large interaction effects among multiple model para- 
meters on the model output (i.e. binary interactions 
terms such as Dm and tertiary interaction terms such as 
Djj k ). This was verified by computing first-order sensitiv- 
ity indices, which account for contribution of each indi- 
vidual parameter (Dj) to the output variance without 



c 

CD 

-t— ' 

2 

Q_ 

D) 
E 

o 
E 
a. 

j/> 

CO 

o 
'c 

CD 
CO 

M— 

o 
c 
o 



c 

CD 

o 
c 
o 

o 



threshold=0.01.LiM 
threshold=0.1uM 
threshold=0.3uM 
threshold=0.5LiM 
threshold=1j.iM 
-|- DMAexp 




Concentration of iAs'" (nmol/mg protein) 
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hepatocytes and the medium) of DMA, while the right panel shows intracellular levels of DMA. Experimental data are from Drobna et al. [53] for 
hepatocytes from a 63 year old white female. 
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taking into account higher-order interactions [57,58]; 
these indices for all parameters were low (< 0.1) at low 
doses (results not shown). At higher doses, k ATG 
(Parameter 1) and (Parameter 6) are the most 

influential. Both these parameters correspond to rate 
constants in the bidirectional reaction of ATG with 
iAs (glutathione conjugation of arsenite and ATG 
hydrolysis). 

A fundamental hypothesis in this modeling formula- 
tion that allows the capturing of the dose-response pro- 
files of arsenic retention and methylation across doses 
(Figure 5), is the introduction of threshold-dependent 
non-linear ("switch-like") mechanisms in the metabolic 
network due to oxidative stress (TR inactivation). This 
assumption is based on findings that signaling motifs 
exhibit biological switches under a narrow range of 
endogenous or exogenous stimuli [59]. This is often 
described by a Hill equation with a large Hill coefficient 
(e.g., kinase cascades [60] and nuclear-receptor pathways 
[61,62]). The large Hill coefficient for inhibition of 
MADG hydrolysis (Table 1 - Parameter 8) points to 
potential "switch-like" behavior of the activation of Nrf2 
due to arsenic-mediated oxidative stress [32]. 

The parametric comparison between human and rat 
hepatocyte TK models for arsenic, presented in Table 2, 
provides insight into factors that affect arsenic metabo- 
lism in hepatocytes. Specifically, the AS3MT activity is 
found to be not a significant factor. Rats have been 
reported to be much faster metabolizers of arsenic than 
humans [51], but based on this study this is attributed 
to other factors. Specifically, in the rat-hepatocyte TK 
model (Figure 1), the protein-bound arsenite (p-iAs) is 
biotransformed to iAs (substrate of methylation 



reaction) at a much higher rate (five orders of magni- 
tude) compared to the biotransformation of iAs 111 to 
ATG (substrate of methylation reaction) of the human- 
hepatocyte TK model (Figure 3). This variability could 
be influenced by polymorphisms related to GSH pro- 
duction in hepatocytes [63] or availability of thiol-con- 
taining proteins to interact with AS3MT [28]. 

This work demonstrates the development of a proto- 
type semi-mechanistic toxicokinetic (TK) model for 
arsenicals in human hepatocytes introducing features 
such as cooperativity and "switch-like" antioxidant 
response. Even though this model is not directly applic- 
able to in vivo systems as a standalone formulation, it 
can be applied to inform macroscopic metabolism- 
related parameters in the PBTK model. On the other 
hand, more experimental studies on arsenicals in human 
hepatocytes will substantially improve model structure 
and can help in characterizing inter-individual variability 
in arsenic metabolism. Currently, the Styblo et al. study 
[51] is the only study in the authors' knowledge that 
reports time course profiles of arsenic methylation in 
human hepatocytes. Furthermore, significant uncertain- 
ties exist in experimental data due to the limitations of 
widely used techniques such as hydride generation- 
atomic absorption spectroscopy (HG-AAS) and high 
performance liquid chromatography-inductively coupled 
plasma-mass spectrometry (HPLC-ICP-MS), where glu- 
tathione conjugated arsenic species ATG and MADG 
have been reported to be degraded to iAs and MMA 
during the speciation analysis in the bile of rats exposed 
to arsenic [64,65]. 

This cellular-level TK model is based on an arsenic 
biotransformation pathway where arsenic-GSH adducts 
(ATG, MADG) are substrates for the respective 
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methylation reactions [27,28]. It should be pointed out 
that arsenic can be efficiently methylated even in the 
absence of GSH [49,66], indicating that arsenic-GSH 
complexes need not be major species in the methylation 
of arsenic. On the other hand, the explicit consideration 
of arsenic-GSH complexes allows the description of a 
hysteresis behavior associated with methylation reactions 
and the stimulating role of GSH in these processes 
(Equations 7-8); so, this mechanism has been selected in 
this model. It should be noted that it is beyond the 
scope of this manuscript to comparatively evaluate the 
different arsenic biotransformation mechanisms. 

Clearly it is, in principle, possible to incorporate into 
this cellular level TK model both the oxidative and 
reductive mechanisms as individual pathways. However, 
currently available experimental data are not adequate 
for estimating the relative contributions of each path- 
way. Development of improved experimental techniques 
for quantifying binding of arsenicals to GSH and thiol- 
containing proteins will allow the estimation of the 
relative contribution of each pathway. Since AS3MT 
coexists in hepatocytes with a number of competing ele- 
ments that affect its action, its activity should be deter- 
mined based on the availability of each of these 
elements. For instance, to study the effectiveness of the 
oxidative mechanism, it is possible to knock out GSH 
biosynthesis in hepatocytes by interfering with the activ- 
ity of Glutamate-Cysteine Ligase (GCL) [67], and expos- 
ing them to various doses of arsenite. Such data sets can 
be used to estimate AS3MT activities along the two 
competing reaction pathways; this type of information is 
necessary in order to extend the mathematical formula- 
tion of the model described here to include both com- 
peting methylation pathways. 

Parameter identification is an important issue in com- 
putational biology since most of the models involve 
more parameters than the available data. The TK model 
was parameterized using data on total arsenic of three 
species (iAs, MMA, DMA) and it was able to capture 
the modes of arsenic retention and methylation in 
human hepatocytes, but was not able to exactly capture 
the time-course profiles from the experimental data. In 
order to reduce the uncertainty associated with this 
issue, sensitivity analysis and testing were conducted as 
a means to identify the relative impact of each para- 
meter on model predictions [68]. Additional time-course 
data (either on intermediate species or under more 
exposure/dose conditions) can improve model perfor- 
mance, and can help obtain additional mechanistic 
insights into the dynamics of arsenicals in hepatocytes. 

Conclusions 

A cellular-level TK model was developed based on a 
recently proposed pathway of arsenic biotransformation. 



This model can describe uptake, retention and clearance 
of arsenicals in human hepatocytes using a semi- 
mechanistic approach. It highlights the key biological 
processes that influence arsenic metabolism by explicitly 
modelling the metabolic network of GSH-adducts for- 
mation [27]. Moreover, comparison of the model struc- 
ture and parameters with a rat-hepatocyte TK model 
[34] highlights the relative roles of different metabolic 
reactions in the methylation of arsenic. Ongoing work 
involves incorporating this cellular-level semi-mechanis- 
tic TK model as a module within a whole-body PBTK 
model of arsenic [44], in order to improve the PBTK 
model parameterization and its predictions [69]. 
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